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Abstract 

The growth of small planetesimals into large planetary embryos occurs much before the dispersal 
of the gas from the protoplanetary disk. The planetesimal - gaseous-disk interactions give rise to 
migration and orbital evolution of the planetesimals/planets. Small planetesimals are dominated by 
aerodynamic gas drag. Large protoplanets, m ~ O.lMm, are dominated by type I migration differential 
torque. There is an additional mass range, m ^ 10^^ — lO^^g of intermediate mass planetesimals 
(IMPs), where gravitational interactions with the disk dominate over aerodynamic gas drag, but 
for which such interactions were typically neglected. Here we model these interactions using the 
gas dynamical friction (GDP) approach, previously used to study the disk-planet interactions at the 
planetary mass range. We find the critical size where GDP dominates over gas drag, and then study 
the implications of GDF on single IMPs. We find that planetesimals with small inclinations rapidly 
become co-planar. Eccentric orbits circularize within a few Myrs, provided the the planetesimal mass 
is large, m > and that the initial eccentricity is low, e < 0.1. Planetesimals of higher masses, 

m ~ 10^^ — inspiral on a time-scale of a few Myrs, leading to an embryonic migration to the 

inner disk. This can lead to an over-abundance of rocky material (in the form of IMPs) in the inner 
protoplanetary disk (< lAU) and induce rapid planetary growth. This can explain the origin of super- 
Earth planets. In addition, GDF damps the velocities of IMPs, thereby cooling the planetesimal disk 
and affecting its collisional evolution through quenching the effects of viscous stirring by the large 
bodies. 


1. INTRODUCTION 

Planets form in protoplanetary disks around young 
stars. Once km sized planetesimals have been formed, 
their evolution is determined by three basic dynami¬ 
cal processes: Viscous stirring, dynamical friction and 
coagulation or disr uption through collisions (see Gol- 
dreich et al. 2004 and references therein for details). 
i'hese dynarnical processes do not include the effects 
from planetesimal-gas interaction in the disk, which can 
be important during the early stages of planet formation 
when gas is abundant (the first f ew Myr, with possib le 
suggestions for longer time-scales, Pfalzner et al.|[20l4 |. 

During the evolutionary phases ot planet formation, 
small dust grains successively grow to large planetary 
embryos. For small size planetesimals, aerodynamic gas 
drag is the dominant effect of the gas. It maintains low 
relative velocities and keeps the orbits circular and co- 
plan ar. Hence, the plan etesimal disk is expected to be 
thin (|Ohtsuki et al.|2002|) . It may also assist in the coag¬ 
ulation and merger ot sm all bodies (]Ormel &: Klahr|20I0' 


Perets & Murray-Glay 20111. Large planetary ernbryos 
[m > U.lAfm) can migrat e due to the interaction with 
the gaseous disk (see e.g. Papaloizou & Terquem 2006 
for a review and references therein.) 

Gas-planetesimal interactions are therefore important 
both for low mass planetesimals and large earth-sized 
or larger embryos. However, there exists an interme¬ 
diate planetesimal mass (IMP) range, of the order of 
TO ^ 10^^ — 10^^g in which gas-planetesimal interactions 
are typically neglected, since such planetesimals are large 
enough to be fully decoupled from the gas, and aerody¬ 
namic gas drag is too weak to affect their gravitational 
dynamics, while they are not sufficiently massive to ex¬ 


ert significant tor que on the ambient gas and change its 
globa l properties (Hourigan & Ward 1984 Tanaka & Ida 
19991. Nevertheless, in this mass range, gravitational m- 


teractions with the disk are not negligible, and the disk- 
planetesimal interactions can be modeled through gas 
dynamical friction (GDF), which dominates over aero¬ 
dynamic gas drag. Here we focus on the effects of GDF 
on IMP and consider its role in their dynamical evolu- 
tion. _ 

Ostriker (1999) showed that dynamical friction in 


gaseous medium exerts a force even at low subsonic ve¬ 
locities. GDF is important for various astrophysical sys¬ 
tems: migration of globular clusters in galactic gaseous 
haloes, in-spiral of binary stars in the companion’s enve- 


Stabler 2010 

Escala et al. 2004b|a Baruteau et al. 12011 ). 

Moreover, it 
comparable t 
km in diamet 
systems with 
tions have be 

was estimated that aerodynamic drag is 
0 GDF for bodies as large as a few hundred 
er. Direct applications of GDF on planetary 
sufhcientlv large eccentricities and inclina- 

en recently studied (ITeyssandier et al. 12013 

Ganto et al.|2013 Rein|2012 Muto et al.|201ip, however 

they deal wit 
IMP range. A 
the literature 
in the non-lii 

h masses ot fully formed planets, above the 
/arious generalizations of GDF are found in 
;; e.g. numerical modeling, circular orbits 

rear regime and others (Sanchez-Salcedo & 

BrandenburgI 20011 Kim & Kim||20U7l |Kim et ai. 1120081 

Kim & Kim 

2009| |Kim 2010p, and more recently |Muto 


ometry, directly applicable to protoplanetary disks. 

Motivated by recent developments both in GDF the¬ 
ory and its application to planetary systems, we explore 
the implications of GDF on single IMP. The implications 
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on binary planetesimals will be discussed in a subsequent 
paper (Grishin & Perets, in prep.). To do so we review 
and com pare the effects of aerod ynamic gas drag (Sec¬ 
tion |2.1[ ) and GDF (Section |2.2[ ) and derive the typical 
size ot planet esim als for which GDF dominates over gas 
drag (Section 2.31. We then calculate the time-scales for 
variation of orbital elements using analytical arguments, 
and study GDF effects numerically (Section]^. Finally, 
we discuss (Section]^ the implications of sum processes 
for the evolution planetesimals. 

2. GAS PLANETESIMAL INTERACTION 

The standard models of evolution of planetesimal 
disks account for various dynamical processes, includ¬ 
ing both physical processes due to planetesimal inter¬ 
actions such as viscous stirring, dynamical friction and 
physical collisions followed by coagulation, as well as gas- 
planetesimal interactions through aerodynamic gas-drag 
on small planetesimals, and planetary migration through 
disk torques acting on large planetary embryos and plan¬ 
ets. 

Although gas-planetesimal interaction has been stud¬ 
ied in the context of small planetesimals, massive plan¬ 
etesimals have been largely ignored as they are decou¬ 
pled from the gas, and their size and velocity disper¬ 
sion is dominated by gravitational interactions. In the 
following, we revisit this argument by carefully exam¬ 
ining processes caused due to the presence of gas and 
discuss their possible implications. We first review the 
basic properties of aerodynamic gas drag and GDF re¬ 
spectively, and then compare both forces and find the 
lower limit of the planetesimal size for which GDF be¬ 
comes dominant compared with aerodynamic gas drag. 
The lower limit is found to be ro ughly in the ord er of a 
few hundred km, compatible with Ostriker| (1999 )’s esti¬ 
mation, and depends on the position in the disk and its 
gas density. 

2.1. Aerodynamic gas drag 

The general drag force imposed on a planetesimal of 
radius R and relative velocity Vrei moving through a 
gaseous medium of density Pg and speed of sound Cg de¬ 
pends on R/\, where A = 1/nF, is the mean free path of 
the gas, n is the number density of the gas, and F is the 
cross section of gas-gas collisions. 

When i? < A, Epstein regime applies where individual 
scattering is considered. The drag force is 


depends on the geometry of the object, but for spher¬ 
ical objects it depends only on Reynolds number, i.e. 
Cd = CniRe). An empirical formula can be u sed for 
CniTZe), fitted fo r the range log^g TZe G [—3, 5] by Brown 
fc Lawyer] ( 2003| ). 


24 

CoiRe) = —(l+0.277^e)°■'‘3-h0.47[l-exp(-0.047^e°■3®)] 
Re 

In the large Reynolds number limit Cu is constant, while 
in the low Reynolds number limit {Re < 1), Cd ~ Re~^. 
From consistency of Eq. 0 and Eq. 0 , the transition 
between Epstein and Stokes drag regimes is i? « (9/4)A, 
and Cd = 24/77.e for Re < 1. 

2.2. Gas dynamical friction 

Gonsider a perturber with mass mp moving on a 
straight line with constant velocity v^ei in a uniform 
gaseous medium with density pg and characteristic sound 
speed Cg. The perturber generates a wake, which in turn 
affects the pert urber. Using linear perturbation theory, 
Ostrikei'l (|1999|) calculated the gravitational drag force 
felt by the perturber. The gas dynamical friction (GDF) 
force is given by 


GDF — 


elAM) 


( 4 ) 


-^rel 


where A4 = Vrei/cg is the Mach number, and 1-{AA) is a 
dimensionless factor given by 


HM) = 


5i"(b4i)--w 

|i»(i-A-) + i»(el) 


M < I 
M > I 

V-rel^ > R m 


Fd = -^FR^PgVthVrel 


( 1 ) 


The force is non-vanishing in the subsonic regime, 
while in the supersonic regime, a minimal radius Rmin 
is introduced to avoid divergence of the gravitational po¬ 
tential (usually taken to be the physical size of the per¬ 
turber, or the accretion radius Cmpjvl^f). The exact 
value of Rmin is not well determined; but it can be fit¬ 
ted through comparison of Eq. with hydrodynam- 
ical simulations, to find a best mting value for Rmin 
( Sanchez-Salcedo &: Brandenburg||I999 1. It is important 
to stress that Ustriker’s original calculation was done 
using a point mass perturber; GDF is a gravitational 
volume force, essentially different from the aerodynamic 
drag, which is a surface force, dependent on the geometry 
of the perturber. 

For small Mach numbers AI <C 1, 


Where Vth = (8/7r)^/^Cs is the mean thermal velocity (for 
Maxwellian distribution). 

For R> X, the gas must be modeled as a fluid. Here, 
the drag force depends also on the Reynolds number 
Re = 2Rvrei/i^m, where Vm = (l/2)ut/iA is the molec¬ 
ular viscosity of the gas. For high Reynolds numbers 
{Re > 800) the gas exerts a ram pressure force, while for 
lower Reynolds numbers Stokes drag is more applicable. 
The drag force is 

Fd = -^CDFPgvj^iVrel (2) 

where Cd is the drag coefficient and Vrei is the unit vec¬ 
tor in the direction of relative velocity. Generally Cd 


I{M) =M^/3 + 0{M^) 

We note that the results of 
qualitatively si milar to 


Ostriker 


( 6 ) 


Kim & Kim (2007) are 


s straight line trajectory 


see fig. (8) of Kim fc kim||2007 ). Gomparison to 


Muto 


et al. (2011 )’s formulae tor slab geometry and validity ot 
both models is discussed in Section |4Tl 


2.3. Comparison between gas dynamical friction and 
aerodynamic gas drag 

Aerodynamic drag is more dominant for small size 
planetesimals and scales as ^ GDF is negligible 

^ For ram pressure regime Co is roughly constant, so this is 
indeed the case. For different drag regimes Co weakly depends on 
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for small sizes, and scales as ^ i?®, since m ~ PmR^, 
where pm is the material density. 

It is therefore clear that there exists a unique value 
Pm,Vrei,'Re,M) for which the aerodynamic drag 
and GDF forces are equal. Moreover, is not depen¬ 
dent on pg. The only independent dimensional parame¬ 
ters are G, pm and Vrei- Dimensional analysis shows that 
i?* scales as i?* oc Vrei{Gpm)~^^'^, where the proportion 
constant depends on the dimensionless numbers TZe, and 
M. Comparing equations ([^ and Q then yields the 
critical size 


i?* = 0.29 


CniRe) 

~I{MT 



'^rel 

^/Gpm 


(7) 


The radial drift of the gas is slightly sub-Keplerian due 
to pressure gradients. For circular orbits, the radial drift 
velocity of the gas is Vgas = I’icVl — p, where 


77 = 


dP/dlna 

2pgvj^ 


( 8 ) 


(Chiang & Youdin 2010| . For flat razor-thin disk, S oc 
, P oc a~'^ we get rj = , where Ho = h/a ^ Cg /vk 

is the aspect ratio of the disk ( Armitag^|2013 ). 

Consider a planetesimal travelling in a circular orbit 
in a gaseous disk. The relative velocity of the headwind 
is Vrei = Vk — Vgas = Kvk- It is Convenient to define 
e = 1 — v^l — 77 « 77/2 -I- 0{Hq). Thus, using ([^, K = e. 
For typical ranges of Hq S [0.01,0.05], the range of e is 
« [10-^4• 10-3]. 

For a concrete example, we consider a disk with simi - 
lar parameters as used by Perets fc Murray-Clay ( 2011[). 
The d isk scaling is adapted from [Chiang fe Coldreich] 
(1997)’s simple flared disk model.^ Let us consider a 
specific power law scaling of disk parameters, with re¬ 
spect to distance, a, to the central star. Disk tem¬ 
perature is T = 120(a/AC/)-3GiL. Molecular weight 
\s p = 2.3mH, where uih = 1-66 • is the hy¬ 

drogen atom mass. This leads to Cg = (ksTj = 
6.63 • 10 ^(a/A[/)- 3 / 34 gj^ . g-i_ xhe disk aspect ratio is 
TLo = h{a)la = 0 . 022 (a/A[/)^/^, where h{a) is the disk 
scale height. 

For full evaluation of the relative velocities we will re¬ 
lax our assumption of circular orbits and let the eccen¬ 
tricity, e, be a free parameter. We show in appendix 
[A| that the ratio of the relative and the Keplerian ve¬ 
locity K = VreilvK alters to iL « + e'^ + 0{eY- 

A more complete treatment of the relative veloci ty be- 
tweeii an ec centric orbit and a planet is given in |Muto| 

’m 7 

CgX is the molec- 


et al. (2011). The Reynolds number is Re = 2RvrIiJv., 
where R is the planetesimal size, p 
ular viscosity, and A is the mean free path. Note that 
alth oug h Eq. 0 appears simple, Cu depends on TZe via 
Eq. (1^, which m turn depends on i?*. 

The top panel of Fig. Hshows the numerical solution of 
Eq. (1^. We see that forlow eccentricities, planetesimals 
with R > 500km are dominated by CDF for most regions 
of the disk. We see that the derivative of the solution 
i?*(a, e) is discontinuous where the Mach number A4 «1. 
It happens approximately where e = Hq, the disk aspect 


R via the Reynolds number. For wide range of Reynolds numbers, 
F£) ~ where 1 < CK(7^e) < 2. 



log [a/AU] 



Figure 1. Numerical solution for equation Q for various eccen¬ 
tricities. Top: Critical size of planetesimal as a mnction semi-major 
axis. Each curve corresponds to different orbital eccentricity e, and 
indicates the critical size R*{a) where both forces are equal. Aero¬ 
dynamic drag dominates for smaller radii, while GDF dominates 
for larger radii. Note log-log scale. Bottom: Estimated Mach num¬ 
ber as a function of semi-major axis for the same orbits. Note the 
logarithmic scale. 


ratio. The origin of the discontinuity is the simplified 
discontinuous behavior of R{A4) near A4 given in Eq. 

€; 

The bottom panel of Fig. shows the dependence of 
the orbit averaged Mach number on a. For Hq <C e, 
Vrei = euif and M. = Vreilcs oc , i.e. as long as 

eccentricity is high, the Mach number is a decreasing 
function of a. For Hq ^ e , Vrei = svr ~ HqVr and 
M. oc so the Mach number is an increasing function 
of a. The transition occurs when Hq is a few times the 
orbital eccentricity. 

In the next sections we will investigate the effects of 
CDF on single IMPs. 


3. GDF EFFECTS ON THE EVOLUTION OF SINGLE 
PLANETESIMALS 


Current studies of gas planetesimal interactions have 
typically bee n restricted to the effects of aerodyna mic gas 
drag forces ( Youdinj 2010 Coldreich et al. 2004). Here, 
we implement the etiects of CDF' for IMPs. 


3.1. Formulation of the problem 

Consider a planetesimal of mass mp traveling with rel¬ 
ative velocity Vrei with instantaneous orbital parameters 
{a,e,I), where a is the semi-major axis, e is the eccen¬ 
tricity and I is the inclination angle. 
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First, we compute the migration time-scale due to 
GDF of a circular orbit, and then compare it to N-body 
simulations with external drag force. For now, we set 
the inclination to zero, thus the problem has two dimen¬ 
sions. Later on we consider inclined orbits and relax this 
assumption. 

Let us consider a typical gaseous disk with Mgas/Mi, ^ 
Hq ^ 0.01 {Mgas/Mi, ~ 0.1 for most massive disks). 
Above this limit self gra vity is non-neglig ible and the disk 


and 


becomes unstable (e.g. [Armitage 2013 and references 


therein). A typical intermediate mass planetesimal has 
a diameter of 200 — 5000 km, and corresponding mass 
of 10^^ — 10 ^® 5 , where the higher end can already be 
considered to be a planeta ry embryo. 


We use Ostriker 


(|1999l)’s linear theory. Deviations 
from linear theory and other models are discussed later 
on. 

The problem of a ccretion of gas by a spherical body 


was first studied by Bondi, H. ( [1952 ). The effective ra¬ 
dius of accretion is Bondi radius tb = 2Gm/c^^. In our 
case, tb/t = S-nGpmr"^/ic^s = 1.25(r/1000fcm)^. Thus 
for masses lower than lO^^g, accretion is negligible. For 
a planetary embryo, we estimate the accretion rate by 
taking 50% accretion efficiency with m ~ 0.57rr^pgCs ~ 
-5“^. It will take Tacc m/m ~ IMyr for a 
planetesimal to accrete its own mass. For simplicity we 
neglect accretion in our simulations. We’ll briefly discuss 
the implications of accretion in the summary. 

3.2. Timescales for the variation of the orbital elements 

Consider a distorting force F = F^r -\- Fg,(p. The vari¬ 
ation of orbital elements can be calc ulated analytically. 
The c hange in the semi-major axis is (Murray & Dermott 
19^ 


da ^ 

mj,— = 2 - , 

^dt y/GMi,{l - e2) 


[FrC sin / -I- ^,^(1 -I- e cos /)] 


and 


( 9 ) 


de 
’ dt 


/a(l-e2) 

GM, 


[Fr sin f + Fg,{cos f + cos E)] (10) 


where a is the semi-major axis e is the orbital eccentric¬ 
ity and / and E are the true and eccentric anomalies re¬ 
spectively. The relationship between the true and eccen¬ 
tric anomalies is tan(//2) = A/(T+~e)7(W-l0 tan(£'/2). 
For small eccentricies, both anomalies coincide. 

The GDF force depends on the relative velocity, hence 
a crucial step in calculating the interaction between gas 
and planetesimals is evaluating the relative velocities. 
The general case is hard to estimate analytically and 
requires avera ging the relative v elocity on orbital pe¬ 
riod (see also Muto et al. 2_011|), but for small eccen- 
tricites (i.e. using equatioiis and the components of 
the force are F). = FQVrei,r and = FQVrei,ip, where 
Fq = AnG^mpPg/Sc^ and Vrei,r and Vrei,!p are given by 
eq’s (Al) and (A2). Neglecting the (1 —term in the 
denominator and assuming p 1, the relative velicities 
are Vrei,r/vK ~ esinf and Vrei,ip/vK ~ —^7/2 — ecos/, 
which leads to a simplied formulae for Eqs. (§ and ([Io| 


a rrir, 


- 6 ^ 008 ( 2 /) - I-t - I) ecos/ (11) 


erur 


[—e{cos{2f)+cos^f)+pcosf] (12) 


averaging over one orbit leads us to (d/a) = —pFo/mp 

and (e/e) = —Fo/2mp, where 27r(A) = Adf . The 
typical time-scale for in-spiral is then 


Anp G'^pgmp 
5.7.106(^^' 


mp Y 
pj V2-10253 j 


yr (13) 


Where po = 3 • 10“'^<? • cm~'^. For m = 2 ■ 

Ta 5.7Myr, which is comparable with the disk lifetime. 
Smaller masses are less affec ted by GDF over the disk 
lifetime (see also Section 13.4| ). Note that is inversely 
proportional to pg. At larger separations, pg decreases 
and GDF is less effective. At lower separations GDF 
stronger, and the rate of in-spiral increases. Thus, Tq is 
an upper limit. 

For small eccentricity, the eccentricity damping 
timescale is much faster than the migration timescale 
(i.e. Te/Ta = 2p = 6iFg ). In addition, the eccentricity 
decays exponentially, i.e. e oc e. The latter is consistent 
with the analysis of eccentrcity waves raised by a planet 


studied by [Tanaka fc Ward 2004 They find that the 


Tfi « H'oTa .For larger ec- 


decay is exponential and that 
centricities, i.e. e > Ffn. the ca lculation is more complex. 
Papaloizou & Larwood (2000) find e cx and 


^/^ 


(14) 


Both the exponential and the power law decay are con- 
firmed in numerical s imulations of type I migration 
(Gresswell et al.||2007|). The connection b etwe en type 


migration and GDh' is discussed in section Af2 
Generally, we expect that the eccentricity to be 
damped faster than semi-major axis. 

The inclination decays rapidly compared with the 
other orbital elements. Orbital inclination changes due 
to the z component of GDF. For an inclined orbit with 
inclination angle /, the normal force is proportional to 
Fz o<i Vz Ivk- For a circular orbit, it is much larger 
than Fg:,, since the ambient gas does not possess any sig¬ 
nificant z component for gas velocity. The inclination de¬ 
cay time is then mdl/dt = ^/ajGM/,Fz cos{ijJ + f). Most 
of the force is applied during the planetesimal passage 
through the bulk of the disk, i.e. where cos{w -|- /) ^ 1. 


/ ^ ^ or 

ri Ta = yv = 2/^ 
la Fz 


2 /"^ = 2 £ 
Ivk 


10 


-3 


(15) 


and we therefore expect that initially inclined orbits 
will be damped on mu ch shorter t ime-scales. For m ^ 

2-10253, Ti ~ 8 -lOV- 

inclined o rbits, rr/r,. 


Rein 


(2012) found that for highly 
I sin (7/2)/sin/ (see Eqs. 10 
and 11 in |Rein||2012| ). In the limit of small inclination, 
I < Hq, it reduces to Tj/va ~ P ^ i/n , which is con¬ 
sistent with Eq. (15) since in section 2.3 we have shown 
that 


£ ~ Hq . We see that generally orbital eccentric¬ 
ity and inclination are damped faster than semi-major 
axis, hence most of the orbits will in-spiral in time-scale 
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comparable to circular orbit given in Eq.(13). For larger 


inclinations the decay time-scale will increase. 

Inclined orbits evolve in three dimensions, hence a ver¬ 
tical structure of the gas density must be introduced. 
We assume a standard Ga ussian vertical structure of the 
form pg ~ exp(—(Youdin 2010), where h is the 
disk scale height. For inclinations much larger than the 
disk scale height, we expect tj to be be larger, since the 
planetesimal spends most of its evolution above or below 
the disk bulk where gas density is low. 

Now that we developed a basic analytic qualitative un¬ 
derstanding of the expected effects of GDF, we continue 
to a more detailed quantitative study of these effects us¬ 
ing numerical simulations. 

3 . 2 . 1 . Numerical set up 

To study the effects of GDF we use an N-body in¬ 
tegrator with a shared but variable time step, using the 
Hermi te 4th order integration scheme following |Hut et ^ 
(]1995 ). In order to include GDF effects we add a fiducial 
GDI' force that mimics Eq. Q. At each step we calcu¬ 
late the additional external acceleration and jerk due to 
GDF. Full description of the calculation can be found in 
the appendix [Bj 


3 . 3 . Results 

In the following we present the results of a single plan¬ 
etesimal evolution and effects of GDF on its evolution 
where various types of orbits and planetesimal masses 
are considered. 


3 . 3 . 1 . Eccentric orbits 

On the top panel of figure we see that all the simu¬ 
lated orbits evolve in a similar manner. This is mostly 
due to the rapid circularization of the orbits; as the or¬ 
bits circularize the relative gas-planetesimal velocities de¬ 
crease, with a corresponding decrease of the Mach num¬ 
ber. 

Eventually, their Mach number decreases to A4 « 1, 
at which stage GDF is highly efficient, leading to a rapid 
loss of the angular momentum. The larger the initial 
eccentricity is, the longer it takes the orbit to reach the 
trans-sonic limit. Generally, the eccentricity damping 
time-scale given by Eq. (14) is compatible with simula¬ 
tions. Actually, Eq. ([l4|) overestimates More rigorous 


derivation could determine Tg more accurately. 

We see that in most cases the planetesimals migrate 
wit hin 3 Myrs, compatible with the estimation of 
in (13). The greater the initial eccentricity of the or¬ 


bit IS, the faster it spirals in. This is due to the en¬ 
hanced gas density {pg oc in the apastron. For 

smaller planetesimals, e.g. of mass ^ Ta increases 

to ~ 200Myr, much more than typical disk lifetimes. 
Note that given the short circularization time of eccentric 
orbits the time-scale for the evolution of the semi-major 
axis is almost independent of the initial eccentricity, and 
orbits with initial eccentricities of e < 0.3 are essentially 
circular through most of their orbital evolution. 

Orbits with initial eccentricity of e = 0.1 circularize 
to ~ 0.02 over 10^ years, and are fully circularized after 
10^ years. For orbits of initial eccentricity e = 0.3, the 
circularization time extend to ~ O.lMyrs. For smaller 
masses, the lowest mass for which the Te is comparable 



Time [Myr] 




Time [Myr] 


Figure 2. Evolution of orbital elements of various orbits starting 
at a = lAU and different eccentricities. The mass of the planetesi¬ 
mal considered here is 2 • Top: Evolution of the semi-major 

axis. Middle: evolution of the eccentricity. Note logarithmic scale. 
Bottom: Evolution of the Mach number. Note logarithmic scale. 

to typical disk lifetimes is ^ lO^^r;, which corresponds to 
500km planetesimals. Less massive planetesimals hardly 
evolve within the disk lifetime. 

3 . 3 . 2 . Mildly inclined orbits 

In figure we show the results for initially inclined or¬ 
bits. In thenrst three panels, orbits with low inclination 
show a fast decay. 

The average ambient gas density experienced by plan¬ 
etesimals at high inclination is low compared with the 
low inclination orbits, hence GDF becomes less effective. 
Highly eccentric orbits, e > 0.3, share the same fate. In 
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Time [Myr] 



Time [Myr] 



Time [Myr] 



Time [Myr] 


Figure 3. The evolution of a 2 ■ planetesimal over IMyr staring on inclined orbits. Solid lines correspond to orbits with initially 

low inclination I = 0.05 rad, and dashed lines represent orbits with initially high inclination I = 0.15 rad. Top left: Evolution of the 
semi-major axis. Top right: Evolution of the eccentricity. Note logarithmic scale. Bottom left: Evolution of the inclination. Bottom right: 
Evolution of the Mach number. Note logarithmic scale. 


special cases (e = 0.3, I = 0.15 rad, dashed cyan line) the 
orbit might gain angular momentum and expand, due to 
tailwind near the apastron that overcomes the headwind 
in periastron. 

For low inclination and circular orbit, the inclination 
drops sharply; a planetesimal on a circular orbit loses half 
its initial inclination already after a few Kyr, comparable 
to the estimated time-scale in Eq. [I^ Orbits with high 
inclinations decay much slower, as expected. 

In Fig. we compare the evolution of inclined or¬ 
bits with those obtained for zero inclination orbits. We 
consider inclinations in the range of / = 0 — 0.15 rad, 
comparable with the disk aspect ratio ~ 2Hq^6Hq, cor¬ 
responding to gas density of Icr and Scr an the apas- 
tron/periastron respectively. The eccentricity is taken 
to be either 0 or 0.1. As can be seen in the first panel, 
the final decay time-scale Tq is the same in all cases re¬ 
gardless of inclination and eccentricity, although for high 
inclination the decay is somewhat slower at first. On the 
second panel, we see that Tg is different in each case. 
There is an order of magnitude difference between orbits 
with initial inclinations of / = 0 and / = 0.05 rad, and 
another order of magnitude difference compared with the 
initial I = 0.15 rad case. Thus, Te is sensitive to inclina¬ 
tion. In the third panel both the circular and eccentric 
orbit decay at the same phase for all inclinations, but the 
circular orbits always have lower inclination since the or¬ 
bit is supersonic on the eccentric orbits and GDF is less 
efficient. Only after the eccentricity becomes negligible 
does Ti retains its higher phase. 

3.3.3. Highly inclined orbits 


In the last section we considered low inclination orbits; 
however, observations of the Solar System and other dy¬ 
namical systems show an ample of evidence for irregular 
orbits, with large and retrograde inclinations. In this 
sections we explore effects on GDF planetesimals in such 
high inclination orbits. 

Gonsider for example circular prograde (/ = 0) and ret¬ 
rograde (/ = 180°) orbits. In each case, the dimension¬ 
less force For prograde orbit, Mprograde ~ 

0.033 and Fpj-Qgj-^d^ = J^progradej'^ — 0.011, while for 
retrograde orbit, J^^retrograde ~ and ~ 

10/90 ~ 0.0012, so Fprograde/Fretrograde ~ 9- ^Ve ex¬ 
pect that the prograde orbit will decay nine times faster 
than retrograde orbit. In Fig. [^we plot the evolution 
of the orbital elements of various initially inclined or¬ 
bits. In the top left panel we see that after 1 Myr, a 
planetesimal on a circular prograde orbit (/ = 0 deg, 
solid black line) has migrated by Auprograde = 0.17AU, 
while the one on a retrograde orbit (/ = 180 deg, blue 
solid line) has migrated by Auretrograde = 0.02AU, hence 
Aciprograde/AdY-etrograde ~ 0.17/0.02 ~ 8.5 , COUSisteut 
with our expectations. 

In the top right panel, only the planetesimal on a co- 
planar circular orbit (dashed black line) was able to cir¬ 
cularize within IMyr of evolution. In the bottom left 
panel we only plot orbits wi th initial inclinations of 20 
deg, since both prograde and retrograde orbits are co- 
planar and do not change their inclination. We see that 
the changes in inclinations are mild. On the other hand, 
inclined orbits with inclination of ~ 9 deg lost their ini¬ 
tial inclinations in less than a Myr. We conclude that 
orbits with inclinations higher than / > 4i/o ~ 12 deg 
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Figure 4. Comparison of the orbital evolution of a 2 • 10^^^ planetesimal for a range of initial conditions. Solid lines correspond to circular 
orbits, e = 0. Dashed lines correspond to initially eccentric orbits e = 0.1. The color index is green -7 = 0, red - low inclination I ~ 2Hq, 
blue - high inclination I ~ 677o. Top left: Evolution of the semi-major axis. Top right: Evolution of the eccentricity. Note logarithmic 
scale. Bottom left: Evolution of the inclination. Note logarithmic scale. Bottom right: Evolution of the Mach number. Note logarithmic 
scale. 


are not significantly affected by GDF. 

3.3.4. Disk density profile 

In the previous sections we considered disks with a 
disk density profile normalized with pg ^ corre¬ 

sponding to a surface density profile of Eg ~ a~°‘ with 
a = 1. In the following we consider other disk profiles. 
In Fig. I^we run simulations of co-planar orbits with ec¬ 
centricity of either 0 (solid line) or 0.1 (dashed line). The 
power law density profiles considered are for a = 1,1.5, 2, 
corresponding to the black, green and blue lines respec¬ 
tively. We have normalized the density profiles at lAU. 

As can be seen in the top left panel, migration occurs 
more rapidly for larger a; the time-scale for migration 
is shortened by factor of ~ 3/2 for a = 1.5 and ^ 2 for 
a = 2. The migration time-scale is insensitive to initial 
eccentricity, since circularization time-scale is rapid. The 
fast circularization is confirmed in the top right panel, 
where the eccentricity decays after few Kyrs. Bottom 
left panel shows the changes in semi-major axis for the 
first lOKyr. We see that after the initial fast decay of 
eccentric orbits, they start to decay in the same manner 
once they have been circularized. Only at later times 
there is difference due to different ambient density closer 
to the star. In the bottom right panel, the changes in 
Mach number are almost identical, consistent with the 
eccentricity decay in the top right panel. In the case of 
fixed normalization at lAU, decreases as the density 
profile is steeper, with virtually no effect to Tg . 

In Fig we run co planar orbits with a = 1,1.5, 2, 
initial eccentricities of e = 0.1, 0.3,0.8. This time, we 
normalize each density profile to po at the periastron. On 


the top left panel we see that the more eccentric orbits 
decay slower since their ambient orbit averaged density 
is lower. Conversely to the previous case, orbits with 
different a migrate and circularize at a different phase. 

3.4. Scaling with planetesimal mass 

Fig.H] shows the dependence of the time-scales for the 
evolution of the orbital elements on the ^ rn~^ mass scal¬ 
ing. We consider time-scales comparable to the gas-disk 
lifetime of up to lOMyr. The slowest evolution time- 
scale, Ta, is in found for semi-major axis; only planetes- 
imals more massive than a few times lO^^p are signif¬ 
icantly affected by GDF. The eccentricity damping is 
more rapid. The lower limit for the mass of planetesi- 
mals which are significantly affected (i.e. the damping 
time-scale of at least one of their orbital elements, Tq, Te 
or Ti, is comparable to the disk lifetime) is m « 8 • 
for e = 0.1 and m « for e = 0.3. Generally Te is 

an increasing function of e, but even for high eccentricity 
(e ~ 0.8), we still get Te < Tq. The inclination time-scale, 
T/, is by far the most rapid. For low inclination I = 0.05 
rad we have the lower limit m « 4 • lO^^p . Generally t/ 
is an increasing function both of initial inclination and 
eccentricity. For J ^ 0.15 t/ is ^ 100 times slower , and 
starting with e = 0.1 adds another order of magnitude 
to Tj. 

It is clear that for single IMPs, GDF damps both incli¬ 
nation and eccentricity efficiently for most mass ranges. 
Planetary embryos of mass > lO^^g also migrate inward 
on time-scales of a few Myrs. 

4. DISCUSSION 
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Figure 5. Comparison between the evolution of planetesimals of m = 2 • with various initial inclinations and eccentricities. Solid 

lines indicate circular orbits, dashed lines indicate eccentric orbits with initial eccentricity of e = 0.1. Black, blue and green lines indicate 
initial inclination of 0,20 and 180 deg respectively. Top left: Evolution of the semi-major axis. Top right: evolution of orbital eccentricity. 
Note logarithmic scale. Bottom left: Evolution of Mach number. Note logarithmic scale. Bottom right: Evolution of inclination. We plot 
only orbits with initial inclinations of 20 deg. The solid cyan line corresponds to an initially circular orbit, and the dashed magenta line 
corresponds to an eccentric orbit with initial eccentricity of e = 0.1. 

Before exploring the implications of GDF for the evolu¬ 
tion of planetesimal disks we discuss the various assump¬ 
tions of which we made use as well as potential caveats 
in the approach taken here. We then also compare the 
effects of GDF to formulation of type I planetary migra¬ 
tion. Finally we discuss the potential role of GDF in the 
evolution of protoplanetary disks. 


non-negligible, where non-linear effe cts tend to decrease 
the GDF force by a factor of a few (Kim||2010). 

Time dependent 3D geometry vs. steady state 
2 D geometry: 


4.1. Validity of assumptions and caveats 

Accretion: _ 

We have shown in section [XT] that gas accretion is ne^ 
ligible for planetesimal radius R < 1000km. In section^ 
we have shown that gas drag is negligible for R < 500km. 
Thus, for intermediate mass planetesimals both gas drag 
and accretion are negl igible. This intermed iate regime 
therefore complements Lee & Stabler (2011)’s analysis, 
where the dominant force considered was due to accre¬ 
tion. For larger radii and fully formed planets accretion 
is the dominant force and one should use models with 


Ostriker (jl999 1’s original derivation made use of time 
dependent perturbation theory, where the perturbation 
is turned on at time t = 0. In a time dependent analysis 
there is a non-zero force in the subsonic regime, how¬ 
ever the surprising result is that this force is time in¬ 
dependent. The reason is that contributions from large 
distances (the far held) are not negligible. However, for 
disks, the emitted sound waves eventually reach the ver¬ 
tical edge h of the disk after time t ^ h/cg = 1/fl ^ lyr. 
Hence, after one orbital period the waves cannot propa¬ 
gate further and the contribution from the far held de¬ 
creases until it becomes negligible. 

In order to tackle the problem 


Muto et al. (20111 


solved the equations of motion in a slab geometry, ihey 


accretion as the dominant drag force (e.g. Lee & Stabler 
20 TTl|Lee & Stahler||2014D. 

Emear regime: 


|Kim (]2010|) has studied the non-linear regime of GDF. 
Tne non-linearity parameter is 


B = 


Grun 


rs 


c2ii(M2 _ 1) 2R{M^ - 1) 


For low Mach numbers it reduces to the condition for 
accretion. Hence we conclude that for planetesimals of 
radius R < 1000km the linear regime considered here 
is applicable. For larger embryos non-linear effects are 


introduced an averaged potential at the scale height of 
the disk, i.e. dt = —Gm^j -|- i /2 _|_ (y/i )2 where 7 
is of order unity factor. They decomposed the solution 
into Fourier components and searched for steady state 
solutions. Moreover, they showed that time dependent 
solutions decay as ^ and therefore after t ^ lyr the 
wake persists in a steady state. 

Nevert heless, such a stea dy state might never be 
reached. Muto et al. (20111 note that the assumption 
for a steady state is only marginally satisfied. Moreover, 
there is uncertainty in the 2D approximation, and the 
force might be altered by a factor of a few. In addition, 
protoplanetary disks tend to be turbulent, which may sig- 
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Figure 6. Comparison between the evolution of planetesimals of m = 2 ■ in protoplanetary disks with with different disk density 

profiles normalized at lAU. Solid lines indicate circular orbits, dashed lines indicate eccentric orbit with initial eccentricity of e = 0.1. 
Black, blue and green lines indicate density profile ~ where a. is 1, 1.5 and 2 respectively. Top left: Evolution of the semi-major 
axis. Top right: evolution of orbital eccentricity. Note logarithmic scale. Bottom left: Evolution of semi-major axis zoomed in on first 
lO'^years. Bottom right: Evolution of Mach number. Note logarithmic scale. 
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Figure 7. Same as Fig. but with normalization at the periastron. Top left: Evolution of the semi-major axis. Top right: evolution 
of orbital eccentricity. Note logarithmic scale. Bottom left: Evolution of semi-major axis zoomed in on first lO^years. Bottom right: 
Evolution of Mach number. Note logarithmic scale. 
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Figure 8. Dependence of the time-scales for orbital evolution ob¬ 
tained from numerical simulations. The shaded region is excluded 
due to gas dispersal after lOMyrs. Solid line indicate semi-major 
axis evolution time-scale, Ta(blue). Dashed lines indicate eccentric¬ 
ity evolution rime scale, Tg, for e = 0.1 (red) and e = 0.3 (green). 
Triangles indicate the inclination evolution time-scale, rj, for a 
circular orbit with initial I = 0.05 rad. note logarithmic scales. 


nificantly affect the perturbation evolution. The turbu- 
lence is parametrized by the turbulent viscosi ty v = acs h 


where a is the Shakura-Sunyaev parameter (Shakura & 
Sunyaev lT973t , not to be confused with disk density scal¬ 


ing power law. Consider a Kolmogorov distribution, i.e. 
the flow consists of self-similar eddies, and energy cas¬ 
cades from the largest eddy to the s mallest one, where 


it is d issipated by molecular viscosity (Chiang & Youdin 
20101. In this case, the largest eddies are of order /q ~ h. 


and the eddy turnover time-scale is tg = Iq/vq = 1/17. 
where Uq ~ Cg is the velocity of the largest eddy. After a 
few to, the gas is well mixed and the density wave starts 
to propagate again from the planetesimal outwards. This 
is equivalent, in some sense, to restarting the problem. 
Since to < I/O, the steady state might not be reached, 
and the pressure wave propagation is restarted every 
fewx to- 

For smaller eddies, the characteristic eddy turnover 
time i^tz ~ (Z/Zo)^^^to- Thus, for a wake with character¬ 
istic length I = Cgt, the relevant turnover time-scale cor¬ 
responding to its length is ti/t = For t <C I/O, 

the perturbation is not affected by the eddy current, and 
only for t ^ I/O the turbulent current of the largest eddy 
destroys the wake, and one therefore needs to consider a 
time-dependent evolution following Ostriker. 

Another possible reason for the steady state not being 
reached is that the planetesimals come back to (nearly) 
the original position after one orbit. For subsonic regime 
the gas has enough time to rearrange itself, but for super¬ 
sonic reg ime subsequ e nt pe rturbations of the disk are 
possible. |Kim & Kim (2007) find the number of interac¬ 
tions with the wake increases monotonically with larger 
Mach numbers. 


Another attractive feature in Ostriker (1999)’s linear 
theory, is that it does not deal with viscosity an d dis- 
sipation, which would affect the model studied by |Muto] 


^^Itisworthnotmgthatthe^roportionaUt^constantforJflagd 

\/aHo and vg = -v/acg 


et al. (2011), since these are second order effects. Our ap- 
proach is therefore complimentary and consider the time 
dependent approach following Ostriker. 

If a is low, and the disk is laminar, we can compare 
betwe en both models. Denoting the GDF forces |Ostriker| 
(|1999|)’s 3D model andlMuto et al.l (|2011|)’s 2D model as 


Fo{M) = FoT{M) 


(16) 


and 


F 2 {M, a) = Fo 


a/8 M<1 
Fo/2 M > 1 


where Fq = AirG'^mpPg/v'^^i (see Eqs. 


respectively, 

(30) and (41) 
timate the ratio between both models 


Muto et al. 2011| for details). We es- 


F2 


M < 1 

2 In + in(l _ 7 ^- 2 ) m > 1 


3a 


For a subsonic perturber, with a ~ 10“^, the dif¬ 
ference is 2-3 orders of magnitude. The situation gets 
better for supersonic perturber, where the second term 
is negligible, and the difference is ~ 2\nVt/rmin- The 
Coulomb logarithm is not well defined, but here InA ~ 
In h(a)/rmin ~ 9 — 10, so the difference becomes smaller, 
only one order of magnitude. We simulated planetesi¬ 
mals of 2 • lO^^g, or with radius of lOOOfcm. The minimal 
radius is the minimum of either the physical radius of 
the planetesimal, the Bondi radius or the non-linearity 
radius. For IMPs Tmin is the physical size the planetesi¬ 
mal, of the order of lOOOkm. The two models are there¬ 
fore marginally compatible in supersonic regime, where 
as a much more significant difference is expected in the 
subsonic regime. 

Shear: 

The derivation for linear regime is valid for homoge¬ 
neous g as. In reality , Kepl erian disks have differential 


shear. Muto et al. (2011) suggest that the shear is 
negligible if the distance from the planet’s semi-major 
axis to the instantaneous co-rotation radius, defined by 
ac = a(l — e)/(l -I- e) is larger the the relevant length 
scale, i.e. |ac — a| ^ h, or e > 2iFo- However, they 
neglected the pressure gradients since they were initially 
interested in supersonic orbits. In our case, for circular 
orbit spiral density wave appears where the relative ve¬ 
locity is sup ersonic and the effe ctive Lindblad resonances 
accumulate (Artymowicz 1993), so the morphology of the 
wave will be distorted on scales comparable to the scale 
height ^ Hqq = h. The overall force will therefore po¬ 
tentially differ by a factor of a few. 


4.2. Connection to type I migration 

In section we have compared between aerodynamic 
gas drag and GDF. As mentioned earlier, these forces 
originate from essentially different physical processes. 
Aerodynamical gas drag originates from difference in the 
pressure due to changes in the flow around the body. It is 
proportional to the area and depends on geometry. GDF, 
on the other hand originates from gravitational interac¬ 
tions between the massive body and the ambient gas, and 
proportional to the mass (squared) of the body. Aerody¬ 
namical gas drag acceleration decreases with increasing 
mass, while GDF acceleration increases with increasing 
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mass. The comparison between gas drag and GDF is 
straightforward. 

In the case of GDF and type I migration, however, the 
comparison and the connections are less trivial. Both 
approaches deal with gravitational interactions between 
the planetesimal and the gaseous disk. 

Up to dimensionless factort]^ the GDF torque is 


is more than two order of magnitudes faster than co oling 


T, 


gdf 


G'^mlpg 


while the type I m igration torque (e.g. Armitage 2013 
Tanaka et al.||2002 ) is 


T, 


migration 


I ^ 


nir 


(17) 


Where Eg = 


2pgh is the gas surface density. Both 
torques depend on the planetesimal mass and disk pa¬ 
rameters on the same way. The ratio between both 
torques is Tcdf/T migration ~ {h/a)~^. The factor h/a 
comes fro m the fact t hat the differential torque scales 
with h/a (Ward 1986). This is because there is an in¬ 
trinsic asymmetry between inner and outer torques. In 
some sense, this is analogous to tidal forces that rise due 
to asymmetry of the gravitational forces induced on a 
solid body. One can better see the relation between these 
approaches by estimating the one sided torque using the 
GDF formulation, where the relevant relative velocity is 
the sound speed and the relevant formula is the super¬ 
sonic one since effective Lindblad resonances resides at 
locations where the relative velocity between t he flow and 
the plane tesimal is equal to the sound speed (Artymow- 


icz 


1993). 


inally, in our study we consider disk-planetesimal in¬ 
teractions on eccentric orbits. Note that the effect of 
eccentricity on planet migration was discusse d in several 
studies such as [Papaloizou fc Larwood 2000[ They find 
that for eccentricity larger than ~ l.l.nn the torque re¬ 
vers es. It has been confirm ed by 
and Bitsch & Kley (2010), but t 


Gresswell et al.l (20071 
re embedded planet is 
always migrating inward. The reason is that the torque 
changes the eccentricity, while the migration rate is cal¬ 
culated from the power, which is always negative for 
isothermal disks. In our models we always see inward 
migration. 

4.3. Implications of gas dynamical friction for planet 
formation and the evolution of protoplanetary disks 

Planetesimal disk evolution: When considering the 
evolution of large and small planetesimals in a disk, the 
evolution of the velocity dispersion is governed by vis¬ 
cous heating of the large bodies by themselves and their 
cooling by dynamical friction, whereas small planetesi¬ 
mals m ostly heat through dyna mical friction by the large 
bodies (|Goldreich et al.||2002|). Introducing gas to the 
system gives rise to additional cooling channels both for 
large and small bodies. For large oligarchs, GDF keeps 
the random velocities low, which prevents oligarch col¬ 
lisions. For small bodies, aerodynamic gas drag is the 
main process which dominates their cooling; such cooling 

^ The dimensionless factors are not necessarily order unity. Some 
of them could be very small, e.g. for the subsonic regime the di¬ 
mensionless factor is proportional to ~ 


of small planetesi mals through inelastic collisions (Gol- 


dreich et al.|2002 ) considered before. Moreover, the drag 


force increases with increasing velocity, so any random 
velocity raised by viscous stirring will be suppressed by 
gas drag. GDF therefore keeps planetesimal disks thin 
and cool much longer than otherwise thought, and can 
not be neglected. 

For the upper tail of large ptotoplanets, the 
mass is close to the isolation mass Miso ~ 


O.O7M0(a/A[/)^(Eg/lOg • where the proto - 

planet has cleared its feeding zone (|Pollack et al. |1996). 
A large protoplanet under GDF force could migrate to a 
new environment where additional-planetesimal swarms 
are available for accretion. Hence the protoplanets do 
not grow in isolation and its growth is not limited by its 
local feeding zone. 

Higher relative velocities can also affect the embryonic 
migration rate and eccentricity and inclination damp¬ 
ing. While initially eccentric orbit is rapidly circularized 
in disks around single stars, planetesimals embedded in 
disks of binary stellar systems can have persistent higher 
relative velocities. The origin of high relative velocity 
can be either eccentric co-planar disk, or inclined orbit 


due t o misaligned stellar companion (Rafikov & Silsbee 


2015). 

Super-Earth / hot Neptune formation: 

The observed distributions of planets in short periods 
suggest that lower mass super-Earths / hot Neptunes are 
more common than gas giants, compared with simple ex¬ 


pectations from population sy nthesis models (Ida & Lin 
2008 Hansen fc Murray|2012 ). In order to reproduce the 


observed mass distribution an over-abundance of rocky 
material is required in the inner (< lAU) region. This 
can be achieved either by radial drift of roc ky material 


in the form of dust or small planetesimals (Hansen & 


Murray! 20121, or by enhanced primordial mihimal mass 


extrasolar nebula (MM EN) surface density profile (Ghi- 
ang fc LaughlinpOld). Other suggestions involve type I 


migration of planets formed in the outer region into the 
inner region, which requires additional physical processes 
to concede with observations (i.e. eccentricity damp- 
ing, tidal friction, pressure maxima trapping, jlnamdar] 
fc Schlichting|[2M4 1. Both suggestions suffer from theo¬ 


retical and observational challenges: migration of small 
planetesimals through aerodynamic gas-drag, might be 
too fast and lead to the accretion of the solid mate¬ 
rial by the star, while enhanced primordial MMEN is 
unlikely and might be unstable (|Inamdar & Schlicht- 
ing 2014). Migrating embryos due to GDF could pro¬ 
vide an additional channel for supply of solid material 
into the inner region in the form of intermediate size 
planetesimal/planetary-embryos. GDF embryonic mi¬ 
gration naturally introduces eccentricity damping, and 
thus alleviates the need for external mechanisms for ec¬ 
centricity damping. Moreover, slightly eccentric orbit 
results in a more efficient migration, hence the migration 
starts at the scale of intermediate size planetesimals. In 
addition, for a configuration of eccentric disk in binary 
star system, higher relative velocity are predominant for 
the entire disk lifetime, and the embryonic migration rate 
might be faster and cause even more favorable conditions 
for generating super-Earths in circumbinary disks. 

These embryos could therefore assist in providing the 


































































12 


required over-abundance of rocky material for in-situ for¬ 
mation of hot Neptunes. In addition, the time-scales in¬ 
volved are comparable with disk lifetime, hence the mi¬ 
gration is not too efficient as to lead to material accretion 
into the star, but its time-scales are sufficiently short as 
to supply the material during the lifetime of the gaseous 
disk. 


5. SUMMARY 

In this study we considered the effect of GDF on single 
IMPs. We find that GDF is the dominant drag force 
affecting the evolution of planetesimals larger than 
for which aerodynamic gas drag effects are negligible. We 
explored the effect of GDF in the linear regime and in 
the mass range where accretion and non linear effects 
are negligible (up to ^ We estimated the typical 

time-scales for the evolution of the orbital parameters of 
the planetesimals due to GDF, and further studied them 
using detailed numerical simulations. The main results 
can be summarized as follows. 

Planetary embryos of mass 10^"^ — 10^^g dissipate their 
inclinations and circularize in less than a Myr, regard¬ 
less of initial parameters. Moreover, they migrate in¬ 
ward on time-scales comparable to the gaseous disk life¬ 
time ^ 5Myr. Such embryonic migration may help 
explain the origin of close-in super-Earth planets, that 
might form and grow in the inner parts of t he pro¬ 
toplanetary disk from such embryos (see also Hansen 


fc Murray 
mals ( m 


2012 


for related issues). Smaller planetesi- 


IW^g) may dissipate low initial eccentrici¬ 
ties/inclinations (eo < 0.3, Iq < O.lrad). Planetesimals 
in the lowest mass range (~ lO^^g) can dissipate low ini¬ 
tial eccentricities (eo < 0.1), and slightly damp small in¬ 
clinations (comparable to the disk scale height). The effi¬ 
cient damping of inclination and eccentricity reduces the 
planetesimal random velocities and assist in keeping the 
planetesimal disk flatter and more circular, exchanging 
the planetesimals kinetic energy into the gaseous disk. 
Such evolution could therefore have implication not only 
for the planetesimal disk structure but also for the long 
term collisional evolution and planetary growth in the 
disk. We conclude that GDF can play an important role 
in the evolution of planetesimal disk, and should be ac¬ 
counted for in the study of the early stages of planet 
formation. 
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APPENDIX 

A. TYPICAL ECCENTRICITY IN THE TRANS-SONIC REGIME 


We wish to estimate the critical eccentricity Cc for trans-sonic regime. The calculation is similar to Muto et al. (|20I1 ) 


where they neglect the pressure gradients of the gas and assume g = 0. The relative velocities ih polar coordinates 








are 
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The implicit assumption in obtaining ( |Al| and (A2 1 is that the large planetesimals are unaffected by GDF on a scale 
of orbital period. More formally, the stopping time is large compared to th e orbital period. Indeed , (Al) and (A2) at 
zero eccentricity are obtained by letting ^17 1 in Eqs. (16) and (17) in Perets & Murray-Clay ( 2 Ull) and taking 

the first order. 

The total relative velocity is 


^rel ■ 


^rel'. 


^rel,tp 


a{e,r]) + b{e,ri) co s f 
1 — 
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Where a{e,r]) = 2 — rj + — 2(1 — ? 7 )^/^(l — and b{e,rf) = 2e — 2e(l — ? 7 )^/^(l — Note that a is 

a dimensionless quantity, not to be confused with the semi-major axis a. In order to acquire the average relative 
velocity, we need to average Eq. ( |A3[ ) on one orbital period: 
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= ^ J Vreiif,e,ri)dt= ^ J v^M,e,r])^df 


Where dt/df = a(l — e^)^/^/(uif (1 -|- ecos/)^). The average velocity is 
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Which can be rewritten as 
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Expanding in powers of e we get 


{Vrel{e))=VK (1 - V^l - 7?)^ + \/1 “ 


0 (e^ 


-, 1/2 


For very low eccentricities e <C e, we get {vreiifi)) = svk- For higher eccentricities e ;g> e we get Vrei 
critical eccentricity Cc for which the flow is trans-sonic is Cc ~ 2 ,Hq = 0.044, consistent with Muto et al. (|2011[). 


(A4) 

(A5) 
evK ■ The 


B. NUMERICAL SET UP 

At each step we first evaluate Cg = 0.022 • and pg = 1.7 • 10“^ at lAU in simulation units. For m ost runs we model 
the gas density with a power-law slope of pg ~ . Steeper slopes are considered in Section 3.3.4 Next, we evaluate 


r] = 19/7(cs/uif)^ and the Cartesian representation of the gas velocity Vx,gas = —'Cr,gas sin/ = —(1 — !a as 

well as Vy^gas = Vr,gas cos f = (1 — ?/)where / is the true anomaly. Finally, we evaluate v^ei = Vp — Vgas , 
the scalar Vrei = I'l’reil? compute the Mach number by Al = u^ei/cg. 

The GDF coefficient Cq d f which is calculated for every case is 47rpgl( Al) where X(Al)is0.51n((l-|-Al)/(l — Al)) — 
Al) in the subsonic regime, and In A -|- 0.51n(l — Al“^) in the supersonic regime. The effective Coulomb logarithm is 
In A = hiCst/Rmin- Taking t = lyr and Rmin to be the planetesimal size, we get a Coulomb logarithm of the order of 
10. In a small interval near Al ^^1 we set T(A1) = 10 to avoid singularity and make GDF function continuous. 

The GDF coefficient is inserted into the acceleration per unit mass da = CcDEVrei/v^ei mass 


* = dt 
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where dI{Ai)/dA4 is Al^/(1 — Al^) in the subsonic regime and —2/(Al(l — Al^)) in the supersonic regime. 
The total acceleration and jerk are a = —iripda and j = —mpdj respectively. 












































